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Abstract:  We  study  a  quasilinear  nonlocal  hyperbolic  initial-boundary  value 
problem  that  models  the  evolution  of  N  size-structured  subpopulations  compet¬ 
ing  for  common  resources.  We  develop  an  implicit  finite  difference  scheme  to 
approximate  the  solution  of  this  model.  The  convergence  of  this  approximation 
to  a  unique  bounded  variation  weak  solution  is  obtained.  The  numerical  results 
for  a  special  case  of  this  model  suggest  that  when  subpopulations  are  closed 
under  reproduction,  one  subpopulation  survives  and  the  others  go  to  extinction. 
Moreover,  in  the  case  of  open  reproduction,  survival  of  more  than  one  population 
is  possible. 

AMS  subject  classification.  92D25,  35A40,  65M06 

1.  Introduction 

In  this  paper,  we  consider  the  following  initial  boundary  value  problem  that  describes 

the  dynamics  of  coupled  size-structured  subpopulations  with  nonlinear  growth,  reproduction 
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and  mortality  rates: 

4  +  (4(:r,P(t))?/).c  +  =  0,  {x,t)  G  (0,a’max]  x  (0,T], 

^  /‘lEmax 

<  gI(0,P{t))uI(0,t)  =  CI{t)  +  '^2  'yI'J/3J{x,P{t))uJ{x,t)clx,  £G(0,T],  (1) 

j=i  ^ 

k  uJ(z,0)  =  uJ,0(a;),  a:  G  [0,smax]. 

Here  uT(x,t),  I  =  1, . . . ,  N,  is  the  density  of  individuals  in  the  J-th  subpopulation  having- 
size  x  at  time  t.  and 

P(t)  =  /  wJ (x)uJ (x,t)  clx 

j= iJo 

is  a  weighted  total  population  at  time  t.  The  function  m1  denotes  the  mortality  rate  of  an 
individual  in  the  J-th  subpopulation,  and  ft1  is  the  reproduction  rate  of  an  individual  in  the 
J-th  subpopulation.  The  constant  parameters  0  <  7/,J  <  1  represents  the  probability  that  an 
individual  of  the  J-th  subpopulation  will  reproduce  an  individual  of  the  J-th  subpopulation. 
The  function  g1  denotes  the  growth  rate  of  an  individual  in  the  J-th  subpopulation,  and 
C1  (t)  represents  the  inflow  rate  of  the  J-th  subpopulation  of  zero-size  individuals  from  an 
external  source. 

The  model  (1)  is  a  generalization  of  several  size-structured  population  models  (often  re¬ 
ferred  to  as  distributed  rate  models)  which  have  been  widely  investigated  in  recent  years 
(see  [8,  9,  15,  16,  18]).  Motivated  by  the  fact  that,  in  addition  to  observable  characteristics 
such  as  size  or  age  of  individuals,  non-observable  genetic  characteristics  may  often  play  a 
critical  role  in  the  development  of  the  individuals,  researchers  in  [8]  presented  the  first  such 
generalization  of  the  classical  Sinko-Streifer  model.  There,  the  population  under  consider¬ 
ation  was  treated  as  being  composed  of  several  subpopulations  with  different  growth  rates, 
i.e.,  there  are  inherent  differences  in  growth  between  the  individuals  of  the  population.  This 
results  in  a  system  of  equations  similar  to  (1)  with  the  parameters  g1 ,  B1  and  rn1  being  in¬ 
dependent  of  the  total  population  (i.e.,  effect  of  competition  is  not  accounted  for).  In  [8]  it 
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was  shown  through  numerical  simulations  that  there  is  a  crucial  difference  in  the  dynamics 
of  the  classical  Sinko-Streifer  models  and  those  of  the  generalized  models.  In  particular, 
the  classical  models  cannot  have  dispersion  of  the  density  of  the  population  in  age  or  size. 
Therefore  the  classical  models  are  in  conflict  with  most  of  the  field  data  collected  by  exper¬ 
imental  biologists  (see  [8]  for  more  details).  In  [9]  an  approximation  method  for  the  inverse 
problem  of  identifying  the  growth  rate  distribution  was  studied  and  convergence  results  were 
presented.  This  method  was  subsequently  applied  [18]  to  a  semilinear  model  where  only 
the  mortality  rate  rnT  depends  on  the  total  population  due  to  competition.  Moreover,  the 
convergence  results  for  the  inverse  problem  were  extended  to  this  setting.  In  [10]  the  inverse 
problem  technique  was  used  to  fit  field  data  (mosquitofish  data  which  attains  dispersion  of 
the  density)  to  the  generalized  linear  model.  The  resulting  data  fit  in  [10]  indicates  that  the 
need  for  such  modification  is  crucial  if  these  models  were  to  be  used  as  prediction  tools. 

When  N  =  1,  problem  (1)  reduces  to  a  classical  nonlinear  Sinko-Streifer  model  that 
describes  the  evolution  of  one  population  with  possible  competition  between  individuals. 
For  the  linear  and  semilinear  forms  of  such  a  model  (where  g  =  g(x)  and  ft  =  ft{x)),  several 
approaches  have  been  developed  in  the  literature  for  establishing  existence-uniqueness  of 
solutions.  For  example,  in  [11,  12,  19]  the  semigroup  of  linear  operators  theoretic  approach 
was  used  to  obtain  such  results.  Monotone  approximations  are  developed  in  [1,  2],  and  upon 
passing  to  the  limit  a  solution  to  the  problem  is  obtained,  whereas  uniqueness  is  obtained 
via  comparison  results.  For  the  quasilinear  case  (where  g  =  g(x,P)  and  ft  =  ft  (x,P)),  the 
well-posedness  has  been  discussed  in  [3,  13],  wherein  completely  different  techniques  were 
used  for  establishing  the  existence  of  a  unique  solution  to  this  model.  In  [13]  the  method 
of  characteristics  together  with  a  fixed  point  argument,  is  employed  to  prove  this  result. 
A  difference  approximation  is  developed  in  [3],  and  upon  passing  to  the  limit  a  solution 
to  the  model  is  obtained.  Then  the  Holmogren  Uniqueness  Theorem  is  used  to  establish 
uniqueness  of  this  solution.  To  our  knowledge,  results  concerning  existence,  uniqueness,  and 
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convergence  of  approximations  for  the  general  quasilinear  case  given  in  (1)  with  arbitrary  N 
are  not  available  in  the  literature. 

In  this  paper,  we  develop  an  implicit  finite  difference  approximation  for  problem  (1). 
Techniques  in  the  spirit  of  those  in  [14,  23]  are  used  to  obtain  existence-uniqueness  of  weak 
solutions  as  well  as  convergence  of  the  difference  approximations.  By  a  weak  solution  to  prob¬ 
lem  (1)  we  mean  a  bounded  and  measurable  function  u(x,  t)  =  ( u 1(x,  £),  u2(x ,  uN (x,  t)) 

satisfying 

P^max  /'H-max 

/  uT(x,  t)(p(x,  t)  clx  —  /  uI'°(x)ip(x,0)  dx 

Jo  Jo 


t  PX  n 


( u 1  p>s  +  gIu 1 ipx  —  m1  u  I(p)  dx  ds 


Jo  Jo 


(2) 


N 


+  /  q?(0,  s)  CT(s)  +  /  ryI,J f3J (x,  P(s))uJ(x,s)  dx  ds 


jo  \  j=i  J o  / 

for  t  G  [0,T],  I  =  1, . . .  ,iY,  and  every  test  function  ip  G  C1((0, a:max)  x  (0,T)). 

The  following  regularity  conditions  will  be  imposed  on  our  model  parameters  throughout 
the  paper:  for  any  /  =  1, . . . ,  N 


(HI)  uI,0(x)  E  BV{ 0,.rmax)  n  Loo(0,a:max)  and  u7’°(a:)  >  0. 

(H2)  mT{x,  P)  is  a  nonnegative  continuously  differentiable  function  with  respect  to  x  and 
P. 


(H3)  Vi  r.  P)  is  a.  nonnegative  continuously  differentiable  function  with  respect  to  x  and  P. 

(H4)  g'(  x,  P )  is  a  twice  continuously  differentiable  function  with  respect  to  x  and  P,  g1  (x,  P )  > 
0  for  x  e  [0,  xmgx),  and  g^x^^P)  =  0. 

(H5)  C1  is  a  nonnegative  continuously  differentiable  function. 

(H6)  sup  /3 T{x,P)  <  U\. 

(x,P)E[0,Xma,yi)  X  [0,Oo) 
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(H7)  For  any  sufficiently  small  8  >  0 


sup 

(sr./’jc-  ()./}  (I  x  ) 

(H8)  w1  is  a  nonnegative  continuously  differentiable  function. 

The  paper  is  organized  as  follows.  In  Section  2,  we  develop  a  numerical  scheme  for 
computing  the  solution  of  (1)  and  prove  the  convergence  of  this  scheme  to  a  bounded  total 
variation  function  satisfying  (2).  In  Section  3,  we  present  numerical  results.  In  Section  4, 
we  show  the  continuity  of  the  weak  solution  under  additional  conditions  on  the  initial  data. 
Concluding  remarks  are  given  in  Section  5. 

2.  Convergence  of  Approximations 

The  techniques  used  in  this  section  are  in  the  spirit  of  those  used  in  [14,  23]  to  obtain 
convergence  of  finite  difference  approximation  to  conservation  laws.  However,  it  is  worth 
pointing  out  that  there  are  some  major  differences  between  problem  (1)  and  a  classical 
system  of  conservation  laws.  In  particular,  the  flux  in  (1)  is  a  nonlocal  nonlinear  function 

/'J-max 

of  the  solution  u1  (i.e.,  g1  =  gI{$,'52j= i  /  wI{x)uIdx )),  whereas  it  is  a  local  nonlinear 
function  in  classical  conservation  laws.  Furthermore,  problem  (1)  is  considered  on  a  bounded 
domain  [0,.Tmax]  with  a  boundary  term  that  is  a  nonlocal  nonlinear  function  of  the  solution 
u.  versus  an  unbounded  domain  R  for  a  classical  conservation  law  system.  In  the  sequel,  we 
shall  show  that  such  differences  result  in  two  problems  that  are  very  different  mathematically. 
In  particular,  it  is  well  known  that  for  a  conservation  law  system  it  is  generally  not  possible 
to  obtain  a  bound  on  the  total  variation  for  the  approximating  solutions,  and  hence  to  obtain 
convergence  one  resorts  to  the  compensated  compactness  method  (see,  e.g.,  [23]  for  details). 
However,  a  bound  for  the  total  variation  of  the  approximating  solutions  of  problem  (1)  is 
established  (see  Lemma  3  in  this  section). 


g  (x  +  6,P)  -  g\x,  P) 


m 


(x,P) 


<  C02. 
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The  following  notation  will  be  used  throughout  this  paper:  Ax  = 


^  and  At  =  — 


my~  =  m1  (xj,  Pk),  wIj=wI(xj)  and  CI,k  =  C/(4)- 


n  m 

denote  the  spatial  and  time  mesh  size,  respectively.  The  mesh  points  are  given  by:  Xj  =  j Ax, 
j  =  0, 1,  2,  •  •  •  ,  n  and  4  =  kAt,  k  =  0, 1,  2,  •  •  •  ,  m.  We  denote  by  u-'k  and  Pk  the  finite 
difference  approximations  of  u1  (xjpk)  and  P(4),  respectively,  and  we  let 

<4  =  »'L, pt)’  4  =  pk)-  ~‘k 

We  define  the  difference  operator 

I,k  I,k 

/  j  ,  \  U-  -  U;  1 

(“>  )  =  1  5  " 

and  the  l 1  and  l°°  norm  of  uI,k  by 

I  a1  ~k\  i  =  \uIjk\  A'' 

Halloo  =  maxJ=o,l,2r.,n  \Uj'k\. 

We  then  discretize  the  partial  differential  equation  in  (1)  using  the  following  implicit  finite 
difference  approximation 

(  I,k+ 1  I,k  I,k  I  .k  I  I,k  I,k+ 1 

1  u;  —  it,-  q-  u-  -  nf  ,  TI  n  ,, 

-A - - J_  +  ^  i-1  +  mJ’V-fc+1  =  0,  1  <  J  <  n 

At  A,t  j  j 

44+1  =  +  ESa  EL  Ai 

l  =  E?=,  EL  «>'4+1  a* 


(3) 


with  the  initial  condition 


1  rjAx 

wj’°  =  7T  uIfi{x)dx,  j  =  1,  •  •  •  ,  n,  I  =  1, . . . ,  Ah 

J  (j—l)Ax 


If  we  define 


=  1  +  ^  <?•’*  +  At 1  <  j  <  n,  I  =  1, . . . ,  N, 

then  (3)  can  be  equivalently  written  as  the  following  system  of  linear  equations  for  uk+1  = 

\vhk+1  Iihk+1  „U+1  ,  2,fe+l  2,/c+l  2,A-+1  ,  W,fc+1  ,  W,fc+1  7.JV,A+11T  ^  pNx(n+ 1) 

L“0  ’  “1  i  i  un  i  “0  ’  U1  '&*  ■  ■  )  un  i  •  •  •  i  u0  i  “1  i  •  j  un  J 

(4) 
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where 


f  =  [Chk  •  Y,Y,  7  VPi'kufk  Ax ,  ui'\ c2'k  \  Y.Y.  7  2'Jti'kufk  Ax , 


J=  1  i=l 


J=1  i=l 


v2'k  v2'k 


,  <?"’*  7 N'J&kuf  Ax.  <’*, 


J=1  i=  1 


and  AA  is  the  following  block  diagonal  matrix 

/  0  0 

0  A2’A  0 

Ak  =  0  0  A3’A 


0  0  0  0  AN'k 


with  the  lower  triangular  matrix 


AT'k  = 


g[;k  0  o 

-ALg1/  d['k  0 

o  ~Wdk  d1/ 


0  0  -^9n- i 


Note  that  using  the  assumptions  on  our  parameters  one  can  easily  show  that  equation 
(4)  has  a  unique  solution  satisfying  uk+1  >  0,  A’  =  0, . . . ,  m.  Next  we  will  show  that  the 
difference  approximation  is  bounded  in  l 1  norm. 

Lemma  1  The  following  estimate  holds: 


£|KA||i  <  (1  + A 


u1'0  i  + 


+  N  wiAf)*"1  IC7’*-1!  At, 


and  thus 


pk  <  p  = 

1  Aa  1  max 


(  N 

max  1 1  w1 1 A  (  (1  +  N  uqA£)TO  V  II u1,0 

=i'~'N  v  h 


+J2  +  N  At)ra“i  |C'/,i_1|  At 


7=1  i=l 
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Proof.  Multiply  equation  (3)  by  Ax,  sum  over  j  =  1,  •  •  •  ,  n  and  /  =  1, . . . ,  N  to  obtain 


£ii«u'+1ii.<£  \\uI,k\\i  +  At  c'-‘-  +  ££  7J’J/5 i'kufk  Ax 

1=1  1=1  L  V  J=  1  i=l 


<  1 1  vJ'k  1 1  \  +  At  CT'k  +  iMoolK’ll 


Y  ii^iii  +  XI  AtcI'k + AtNYl  ii/5Jii0oii“J,i1 


<5^11  Ur'k  1 1 !  +  At  Y  C''k  +  AtN  J ™\r  1 1  /5J 1 1  CO  Y 


Since  max  /id (a:.  P)  <  oq,  it  follows  that 


X>J’A+1||1< 


<  (1  +  Nu>iAt)  Y  ll^lli  +  At  £|c"l, 


which  implies  the  estimate.  □ 

We  then  establish  an  l°°  bound  on  the  difference  approximation. 

Lemma  2  Assume  that  At  is  chosen  to  satisfy  oqAt  <  1.  Then  we  have  the  estimate 


«/,A:||oo  <  max 


1  —  uizAt 


IIC'IU  +  ^eL  ll«,'‘-1ll. 


where  a.\  <  g1  { 0,  P),  I  =  1, . . . ,  N. 

Proof.  We  first  note  that  if  max  ?/[J:+ 1  occurs  at  the  left  boundary,  then  from  the  second 


equation  of  (3) 


9ok\u o’fc+1|  <  \Chi\  Put 


Otherwise,  suppose  that  for  some  1  <  j  <  n,  uf  +  =  maxw7,  +  .  Then  from  the  difference 

J  i 

equation  (3)  we  have  that 


/i  .  aj  i ,k  ,  At  Ik  Ik+1  At  Ik  Ik+1  Ik 
' 1  1  A/ 1,1  j  1  T7  fli  )ui  TT-Vi  ",  • 
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Rearranging  terms  and  using  the  inequality  <  «^A:+I ,  we  find 

J,k  _  I,k 

m  i  a  <  I,k\  I,k+1  ,  a  ,  @j— 1  J,fc+1  ^  /,/: 

(1  +  Atm •  )u-  +  At  — — : — - — «_•’  < 


Arc 


—  a 


Hence,  by  (H7)  we  obtain 


(1  —  lu2A t)uTfk+1  <  U  :'k  <  max  u\'k , 


which  implies  the  estimate.  □ 

Multiplying  equation  (3)  by  wj,  summing  over  j  =  1 , . . . ,  n,  /  =  1,...,1V,  and  using 
Lemmas  1-2  one  can  easily  show  that  there  exists  a  c  >  0  such  that 

I  pk+ 1  pk 


At 


<  c. 


(5) 


The  bound  (5)  will  be  used  in  the  proof  of  the  next  lemma  where  we  show  that  our  approx¬ 
imations  Uj'k  have  bounded  total  variation.  This  result  plays  a  crucial  role  in  establishing 
the  subsequential  convergence  of  the  difference  approximation  (3)  to  a  weak  solution  of  (1). 
We  remark  again  that  such  a  bound  is  not  possible,  in  general,  for  a  system  of  conservation 

laws  (see  [23]). 


Lemma  3  Assume  At  satisfies  max{wi,w2}  At  <  1.  Then  there  exists  a  constant  c  = 
c(max  ||mj,0||bv,  max  ||C,/||ci(o,t))  such  that  for  all  k  =  1,  •  •  •  ,  m,  \\Df  (■ uI,k )  ||i  <  c,  /  = 

l,...,iV. 

Proof.  Set  r)j,k  =  Df  ^  Uj'k^j  and  apply  the  operator  Dh  to  equation  (3)  to  get 


I,k+1  I,k  /  I,k  /,fc+l  I,k  I,k+ 1- 

<h  -%  ,D-(V,  «t  |+O-(myV^)  =  0.  2  <  j  <  n 


At 


Ax 


and  for  j  =  1  we  have  that 


J,fc+1  I,k  i  /  /,A-+1  I,k+ 1  J,fc  /,A' 
?71’  —  1  /  tq  —  u  0  !  —  ’ 


ui  -“o 


At  At  \  A  a; 

1 


Ax 


Aa; 


At 
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Multiplying  each  equation  by  Axsgn (77b  +  ),  using  the  fact  that  — ry  ,/csgn(?/j’  +  )  >  —  [r)I,k\, 
and  summing  over  the  indices,  j  =  1,  2,  •  •  •  ,  n,  we  find 


jl,k+l  i  _ 

I  yjl 5^  1 1  ^ 


At 


+E 

j=i 


7,fc  7,/c+l  _  7,7  ,7,7+1  \ 


D.  j  9j  +  j 


,7,fc  ,7,7+1  \ 


sgn(??A’A+1)  Arc  <  0, 


where  we  set  mT(j'k  =  0  and 


D - 


/  7,7  7,7+l\ 

{%  u0  ) 


7,7+1  7,7 

-u0  —  U, 


0 


At 


Now,  simple  calculations  yield 


Ea, 

3= 1 


,  7,7  7,7+1  7,7  7,7+1  \ 

+  -  gj-x^-x  )sgn{J,M]Ax 


Ax 


J 


Thus, 


> 


Ea- 

3=2 


-  £1- 


j-i  7,7+1 
Ax  “J-1 


sgn(//J‘/'  !  Ax  + 


7,7+1  7,7 

«o  —  «, 


0 


■sgn(?h’A+i) 


At 


+D : 


(9,u)  «;'‘+1Sgn(17;-t+'). 


^7,7+1  i  _ 

I  y-jl 1 1  ^ 

^  7  \ 

<  max(D/{  yjr,J  j 


m „• 


+  max,  l-D/j  {Dh  (gj’k) 


mj’A)  |  ||«/’fe+1 


7,7+1  7,7 

«o  -  «o 


At 


£>: 


(^) 


fit 


7,7+1 

0 


From  Lemmas  1-2,  it  suffices  to  obtain  a  bound  for  the  term 
consider  the  boundary  condition 


N  n 

£/oku ok+1  =  CT'k  j  r'J  ufk  Ax. 

J=  1  j=i 


7,7+1  7,7 

"*0 


At 


(6) 


To  this  end, 


10 


Then,  using  equation  (3)  together  with  summation  by  parts  we  obtain 
9o 


A*  -  A*-'  \  c*  -  c*-' 


At 


u  o  - 


J,k  J.k—  1 

"  i  "j 

At 


At 

■  ijk  -  if  r 
At 


J,k- 1 
UJ 


N  n 


<EE  -Pfk  (Dh  (af  'uf)  +  rnfk  +/3Jp(x,P) 


■i  i  /-  i 

N  n  r 

sEE 

j=i  j=i  L 
N 


( D~{0fk)gf  1  -  fifmf  x)  uf  +  /^(x,  P) 


Ax 

pk  _  pk—1 

At 

pk  _  pk—1 


J,k- 1 


Ax 


At 


u 


J,k- 1 


Ax 


,  \  A  oJ./c  J,k— 1  J,/c 
i  E  V  A,  »r,  • 
j=i 

where  P  is  between  P*_1and  Pk.  Hence,  using  the  bound  given  in  (5)  we  find 
9o 


At 


7,/c+i  _  I,fc\  (  I,k  _  I,k—1  \  r<l,k  _  ryl,k- 1 

j,k  i  uo  uo  \  p  I  i  ,.k  ^ 


At 


Uq  - 


N 


J= 1 


<  E  SUP  D~h  (  f  )„f  !  +  sup  |;3f  m,f -1 


At 


iV 


E<  sup  /j)  n«Jvfc  xi 

J=1  £[0,2?max]  X  [0, -Pmax] 


AT 


+  e  i  Pik i  I  c j’^_i  |  +  e  su 


j= i 


N 

\ 

sup 

L=1  J 


/3 


L,k—1 


uL'k~l  u 


Now,  Lemmas  1-2  imply  that  there  exists  a  constant  a  >  0  such  that 


I,k+1  l.k 


Uq  ~  U 


0 


At 


<  a. 


Applying  this  bound  to  (6),  we  conclude  that  there  exists  a  constant  ujq  >  0  such  that 


,/,A-+l  1 1  1 1  ~~.Uk 

1  1,1  <^ll»£+1lli  +  ^, 


At 


and  the  result  is  established.  □ 


The  next  result  shows  that  the  difference  approximations  satisfy  a  Lipschitz-type  condi¬ 
tion  in  t. 
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Lemma  4  Assume  At  satisfies  maxjwi,^}  At  <  1.  Then  there  exists  an  A  >  0  such  that 
for  any  m  >  p 


n  I,m  I,p 

",  ", 


E 

3= 1 


At 


|  Aar  <  A(m  —  p),  I  =  1, . . . ,  N. 


Proof.  Summing  the  first  equation  in  (3)  over  j  and  multiplying  by  Aar  we  obtain 


n  7,/c+l  I,k 

"j  ~  "j 


E 

3= 1 


At 


I  Aar  = 

3= 1 


7,/c+l  7,/c+l  I,k  I,k 

Pk  ",  "j-1  _  1.1.  ■  i  9j  9j- 1  _  /./.  /./.  ■  ! 

s,-i  Ax  ",  Ax  mi  ", 


<  sup 
j 


I,k  I,k 


n  7,A:+ 1  7,/c+l 

tA  a _ ~  J-1  |  Aar  +  sup  rj  ~  9j~l  -l  11,,^+r 


i=i 


Aar 


Aar 


mj  Ml®  ’  111 


<  suP(i,p)G[o,imax]x[o,pmax]  \9T{x,P)\  ||?+'+1||i  +  u2  |  «M'n|  i  <  A. 


Hence 


n  I,m  I,p 

",  -  V 


E 

i=l 


At 


A-r<  EE 


m  n  /,/c+l  „  J,/c 

J  J  -|  Aar  <  A(m  —  /?).□ 


k=p  j= 1 


At 


Aar 


Following  [23]  we  define  a  family  of  functions  {4+^}  by 
uL,At(x’  t)  =  uIj'k  for  ®  G  A?  I  e  [4-1,4),  j  =  !,•■■,  n,  A;  =  1, m,  J  =  1, . . . ,  iV. 

Then,  the  set  of  functions  {t/+ At}  is  compact  in  the  topology  of  L1((0,  armax)  x  (0,  T)),  and 
we  have  the  following  lemma. 

Lemma  5  For  I  =  1, . . . ,  N  there  exists  a  sequence  {U^x.  At. }  C  {t+xA/}  which  converges 
to  a  -BV+0,  armax]  x  [0,  T])  function  uT(x,t)  in  the  sense  that  for  all  t  >  0 

^max 

/  I Uixi, At,  (x,t)~  A ( X ,  t)  I  rfx  ->■  0, 


and 


f*T  PXmax 

Jo  Jo 


1  .'vrj.Mj  {x,  t )  -  n/(ar,  t)  |  clx  dt  ->•  0, 
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as  i  — >  oo.  Furthermore,  the  limit  function  satisfies 


u 


I  |.BV([0,; 


i^max]  X  [0,T]) 


<c(\u 


/,0| 


BV, 


^Wc^om))- 


Proof.  The  result  follows  from  Lemmas  1-4  and  the  proof  of  Lemma  16.7  (page  276)  in  [23]. 

□ 


The  next  theorem  will  show  that  the  limit  function,  u  =  (u1,  ir. ....  uN),  constructed  via 
our  difference  scheme  is  actually  a  weak  solution  of  problem  (1). 

Theorem  6  Any  limit  u(x ,  t.)  =  ( u 1(x,  t ),  u2(x ,  t), . . . ,  uN ( x ,  t ))  defined  in  Lemma  5  is  a  weak 
solution  of  (1)  and  satisfies 


N 


P{t)  <  max 

1=1 . N 


1=1 


<  max  1 1  it/ 1 


Dwi  NT  ' 


N  N  rT  N 

Y  ii^°iii+E  /  ^{T-s)ci{s)ds 

i=i  i=i  / 


=  p 


and 


0^2 T  II  /,0| 


I  "  I|l-((0!^x(0,t))  <  max  ^  e  2  \\u 
where  a  \  <  g1  fid,  P)  for  P  G  [0,  P],I=1,. . .  ,N. 


Halloo  +  ^2I=1  ||«/(f)||i 


Cn\ 


Proof.  This  result  can  be  easily  established  by  using  similar  techniques  as  in  the  proof  of 
Lemma  16.10  (page  279)  in  [23].  □ 


The  following  theorem  guarantees  the  continuous  dependence  of  the  solution  { uJfk }  to 
(3)  with  respect  to  the  initial  condition  u1,0. 

Theorem  7  Let  {ufk}  and  {ufk}  be  the  solutions  of  (3)  corresponding  to  the  initial  condi¬ 
tions  uj’°  and  ulf°,  respectively.  Then  there  exists  a 

o  =  (max  I  If/’*  lb,  max  llu^Hoo,  max  1 1  Hr  (uI,k)  lb) 

k,i  k,i  k,i  n  v  J 
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such  that 


TV 


TV 


Y,  ||«/,A'+1  -  u^+%  <  [1  +  (wi  +  amax|y  ||)At]  || u1*  -  u^k\\i 


1=1 

for  all  A:  >  0. 


/=i 


Proof.  Let  vi'  =  u,'  —  u \'k  for  0  <  k  <  m  and  0  <  j  <  n.  Then  v\'k  satisfies  the  following 


(  /,/c+l  I,k 

1  ~  j  -  +  D-k  (g’%^  -  sf  «tl+1) 


At 


.  I,k  I,k+ 1 

+mj  Vj 


/  I  .k  ^  I .k\  ^/,/c+l  /~v  -I  *  / 

(rrij  —  rrij  )u-  =  0,  1  <  j  <  n 


(7) 


I.k  I,k-\-l  ~I,k~I,k-\- 1  \~^n  T  t  nJ,k  J,k  a 

'  ^  =  Ej=i  Ei=i  r  ’  A 


fl'o  "o 


do  "o 


+  EL.  el.  iw  («“  -  A")  A*, 


where  <L’  =  gI(xj,Pk),  and  similar  notation  is  used  for  the  rest  of  the  parameters.  Multi¬ 


plying  each  equation  of  (7)  by  Ax  sgn(v[’*+1),  summing  over  j  =  1, . . . ,  n,  /  =  1, . . . ,  N,  and 


using  the  following  fact: 


En  r~\ —  (  I.k  I.k-\- 1  ^I.k~I.k-\-l 

j  i  °h  [H.i  uj  ~  9j  ", 


j  sgn(rAfc+1)  Ax  >  -go"\v, 


I,k  |  I,k+1 1 


0 


+el.at 1  («L  -  T*)“7+1  sgn^y+'iA*, 


we  get  that 

TV 

E 

i=i 


V 


I,k+1  ^  ^ 


TV  n 


At 


<  -  E  E  at  [<4*  -  )“f +1J 

i=i  j  i 


tv 


TV  n 


+  S#o’ 


A 


0 


1=1 


;  -  //A-4 ■  'sgut/j-1  ;A.r 

/  I  i-  I 


TV  n 


—  'mj’fc|"j,fc+1|  A;r- 

/=!  7=1 

Now,  using  assumption  (H4)  we  can  easily  obtain  the  following 


-  EE  el.  at  I  (st*  - 


,1.*  M’MMd+l 


sgn(A,A:+1)A.T 


<  (ci  max  ||'U/,A:+1||i  +  c2  max  ||fi/,A:+1||oo  +  c3  max  || Dh  (uI,k+1 
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Furthermore,  it  follows  from  the  second  equation  of  (3)  and  assumptions  (H2),  (H3)  and 
(H6)  that 


E/= 1 9o 


hki^hk+li 


V, 


~  E?=i  E”=i«*  - 


o 


J3- 


I,k 


I,k+1\ 


-  E?=i  E”=i  <  Wl  Ef=i  ll^lli 

+  (C4  max  1 1  u I,k+1 1 1 1  +  C5  max  1 1  fi/,A+1 1 1  ^ )  I  pk  —  pk 

i  i 

Hence,  choosing  a  >  0  so  that 


(ci  +  C4)  max  ||uA+1||i  +  (C2  +  C5)  max  ||h/c+1||i  +  C3  max  \\Dh  (fiA+1)  ||i  <  cr, 


we  obtain 

N  N 

||u/,A+1||i  <  [1  +  (uq  +  a  max  |  \wT\  |oc)Ai]  ||u/,A:||i, 

1=1  1=1 

which  implies  the  theorem.  □ 

Next,  we  prove  that  the  BV  solution  defined  in  Lemma  5  and  Theorem  6  is  unique.  To 
this  end,  assume  that  P(t)  G  C1(0,T)  and  B1  (t)  G  (7(0,  T)  are  given  functions  and  consider 
the  following  initial-boundary  value  problem: 

u\  +  P(t))u'),-  +  m1  (x,  P(t))uJ  =  0,  {x,t)  G  (0, rcmax]  x  (0,T], 

<  gI(01P(t))uI(01t)  =  BI(t)t  t  G  (0, T],  (8) 

k  U^x,  0)  =  uIfi{x),  x  G  [0,  xmax]. 

Then  one  can  easily  show  that  (8)  has  a  unique  weak  solution  (note  that  this  system  is 
uncoupled  and  has  a  local  boundary  condition).  In  fact,  a  weak  solution  can  be  defined  as 
a  limit  of  the  finite  difference  approximation  (3)  with  the  given  numbers  Pk  =  P(tk )  and 
uniqueness  can  be  established  by  using  a  similar  technique  as  in  [23,  Page  282],  Hence,  the 
finite  difference  solution  to  (3)  with  given  numbers  Pk  =  P(4)  and  BI,k  =  B1  (tk)  converges 
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to  the  unique  solution  of  (1)  with  the  given  P  G  C'1(0,T)  and  B1  G  (7(0, T).  In  addition, 
from  the  proof  of  Theorem  7  we  can  easily  show  that  if  {uk}  and  {uA’}  are  the  solutions  to 
(3)  corresponding  to  given  functions  (Pk,BI,k)  and  ( Pk\BI,k ),  respectively,  then  we  have 


E  \\vI" 


<  >  \W 


+  a  At  |  Pk  -  Pk  |  +  At^2\BT'k  -  BT'k  |, 


where  vLk  =  uT'k  —  uI,k.  Equivalently, 


vI,k  i  < 


^lK’°lli  +  ^  a\Pl  -  Pl\  +Y^\Bl'1  ~  B1'1]  At. 


Now,  since  from  Theorem  6  for  I  =  l,...,iV,  {UAxAt}  converges  to  upx.t.)  and  {UAxAl} 

converges  to  upx.t)  strongly  in  C([0,  T];  L1(0, 2;max)),  taking  the  limit  in  (9)  we  obtain 

N  N  rt  r  N 

^||«j(f)||i  <  ^||'«/(0)||i+  /  a\P(s)-P(s)\  +Y/\BI(s)-BI(s)\  ds ,  (10) 

7=1  7=1  ^  L  7=1 

where  u(x,t ),  u(x,t )  are  the  unique  solutions  to  (8)  given  (P(£),  B1  (t))  and  (P(£),  B1  (£)), 
respectively,  and  vpt )  =  up-.t.)  —  up-,t).  Then,  applying  the  estimate  given  in  (10)  for  the 
corresponding  solutions  to  (8)  where 

^  pX  max 

P(t)  =  /  wI(x)uI(x,  t)  dx, 


7=1  J0 


J= 1  -70 


/*max 

PJ(f)  =  C,/(f)  +  /  7/,j/5 J (x,  P(t))uJ (x,t)  dx, 


7=1  Jo 


/*max 

P(f)  =  /  w1  (x)  upx,  t.)  dx , 


j=i  -70 


/*max 

Bpt)  =  Cpt) +  ^2  /  7 I,J f3J (x,  P(t))iiJ (x,t)  dx 


are  defined  in  Theorem  6,  we  obtain  the  following  result. 

Theorem  8  Suppose  that  u  and  u  are  two  bounded  variation  weak  solutions  of  (1)  corre¬ 
sponding  to  initial  conditions  u°  and  u°,  respectively.  Then 


i  v  1  v 

^2  \\uIP)  -  «7(^)lll  <  e[(<T+WmaX/  Halloo)  max,  IhhlooT^1]*  ||^(0)  _  up 0)11!, 
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which  means  that  the  bounded  variation  weak  solution  to  (1)  is  unique. 

Hence,  from  Theorem  8  it  follows  that  the  finite  difference  solution  converges  to  the 
unique  bounded  variation  solution  of  (1). 

3.  Numerical  Results 


In  this  section,  we  provide  some  numerical  results  that  corroborate  the  convergence  theory 
presented  in  Section  2  and  demonstrate  the  effect  of  the  probability  function  Y'J  on  the 
dynamics  of  this  system.  For  the  rest  of  this  section  we  assume  that  amax  =  1,  N  =  2,  and 
the  weight  functions  w1  =  1.  This  implies  that  P  =  [u1(a:,t)  +  u2(x,  t)]  dx.  We  choose 

the  parameters  g1,  ft1 ,  m1,  and  C1  as  follows: 


gT{x,P)  =  g^k1  f1  {P){\  -  a),  ,/?J(.'r,P)  =  /3q  ( 1  -  kI)fI{P)x ,  mT{x,P)  =  mJ(P),  C1  (t)  =  0, 


where  k1  G  (0, 1),  /  =  1,  2.  Hence,  our  model  reduces  to 

u\  +  g^k.1  f1  (P)((l  -  x)u)x  +  m1  (P)uI=  0,  x  G  (0, 1],  t>  0, 

<  glk1  u1  (Q,t)  =  ~  kJ)  J  xuJ(x,t.)clx ,  t  >  0,  (11) 

k  u^x,  0)  =  uI,0(x ),  x  G  [0,1]*  I  =  1,  2. 

Integrating  (11)  and  multiplying  (11)  by  x  and  integrating  once  again,  we  readily  obtain  the 
following  system  of  differential  equations: 

f  M'  =  EL^/?“(1“*J)/J(P)Qj-ra'(F)P'  I  =  hX  (12) 

1  (Ql)  =glk’J'(P)(P' -Q’)-m’(P)Q'.  1  =  1.2, 

where  P1  =  uJ(x,  t)dx  and  Q1  =  xuJ(x,  t)dx ,  /  =  1,  2.  In  our  numerical  simulations  we 
have  used  fl{P)  =  e~°-5P,  /2(P)  =  e~0-lp ,  ml{P )  =  0.3 P,  m2(P)  =  2P/(1  +  P),  k1  =  0.5, 
k2  =  0.7,  /3q  =  1,  S'q  =  1,  t  G  [0, 10]  and 


«1,0(rr) 


5  a:  G  [0,  0.3] 
0  a:  G  (0.3, 1] 


u 


2,0 


a:) 


8  a’  G  [0,  0.2] 
0  x  G  (0.2,1] 


To  test  our  code  we  compute  (P2,  Q1),  I  =  1,2,  the  solution  of  (12)  using  a  4  —  5,/l  order 


Runge  Kutta  routine  available  in  MATLAB  5.3.  We  then  compare  the  total  population, 
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(PJ)A‘C  (t),  and  the  total  biomass,  ( QT)Ax  (£),  that  result  from  the  finite  difference  system 
described  in  Section  2  (with  Ax  =  0.02,  and  At  =  0.01)  to  the  numerical  solution  (P1 ,  Q1^  of 
the  differential  equation  system  (12)  that  results  from  the  Runge  Ivutta  routine.  In  Figures 
1-2  we  present  the  differences  (PT)Ax  (f)  —  PJ(f)  and  (Qj)A'c  (f)  —  QJ(f),  respectively.  The 
figures  indicate  that  the  finite  difference  approximation  provides  a  good  approximation  to 
the  solution  of  (11).  Next  we  present  two  different  dynamics  for  the  total  population  and 
total  mass  depending  on  the  choice  of  Y'J . 

3.1.  Closed  reproduction  case 

In  this  case,  we  assume  that  reproduction  is  closed  under  subpopulations,  that  is,  individ¬ 
uals  in  the  J-th  subpopulation  only  produce  individuals  in  the  J-th  subpopulation.  Hence, 
7/,J  =  1,  /  =  1,2,  and  xI,J  =  0,  /  ^  J.  Solving  (11)  we  present  the  total  populations 
(Pr) Ax  (f),  I  =  1,2,  and  total  mass  (Qt)Ax  (f),  /  =  1,2,  in  Figures  3-4,  respectively.  Note 
that  in  this  case  the  second  subpopulation  goes  to  extinction.  Similar  phenomena  have 
been  known  to  occur  in  different  structured  and  non-structured  population  models,  where 
the  surviving  subpopulation  is  often  referred  to  as  the  fittest  among  the  others  (e.g.,  see 
[4,  17]). 

3.2.  Open  reproduction  case 

In  this  case,  we  assume  that  reproduction  is  open  under  subpopulations,  that  is,  a  sub- 
population  of  type  /  also  reproduces  individuals  of  type  J.  For  our  numerical  example  we  set 
7/,J  =  |,  /,  J  =  1,2.  However,  we  have  performed  many  other  numerical  experiments  with 
different  positive  7 /,J,  and  the  dynamics  are  essentially  the  same.  In  Figures  5-6  we  present 
the  total  subpopulations  (^pT^Ax  (t),  /  =  1,2,  and  total  mass  (Qt)Ax  (t),  respectively.  Note 
that  in  this  case  both  populations  survive  and  approach  an  equilibrium. 
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Figure  1:  The  difference  between  the  total  population  resulting  from  the  finite  difference 
scheme  and  the  total  population  resulting  from  solving  the  differential  equations  system 
using  Runge-Kutta  routine. 


Figure  2:  The  difference  between  the  total  mass  resulting  from  the  finite  difference  scheme 
and  the  total  mass  resulting  from  solving  the  differential  equations  system  using  Runge-Kutta 
routine. 
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population  total  mass  CO  population  total  number 


1.2 


Population  #1 
Population  #  2 


The  computed  total  population  (PJ)A‘C  (£),  /  =  1,2,  on  the  interval  [0,30] 


t 


Figure  4:  The  computed  total  mass  (Qt)Ax  (£),  /  =  1,2,  on  the  interval  [0,30]. 


population  total  mass  population  tote 


Population  #  1 
Population  #  2 


:  The  computed  total  population  (PI)Ax(t)>  I  =  1,2,  on  the  interval  [0,30] 


t 


Figure  6:  The  computed  total  mass  (Qt)Ax  (£),  /  =  1,2,  on  the  interval  [0,30]. 


4.  Regularity  of  Weak  Solutions 


The  goal  of  this  section  is  to  show  that  the  bounded  variation  weak  solution  of  problem  (1) 
is  continuous  provided  that  the  initial  distributions  u 7,°,  I  =  1, . . . ,  N,  satisfy  the  following 
compatibility  conditions: 

N 

(H9)  Let  P°  =  ^2  /0Xmax  w1  (x)  uI'°(x)dx,  gI'°(x )  =  g^x.P0),  f3I,0(x)  =  j3 1  {x,PQ)  and 
1=0 

m1'0  (x)  =  m1  (x,P°),  I  =  1 ,N.  Assume  that  the  function  u1,0  is  a  nonnegative 
continuous  function  and  satisfies 

^  pXmSix 

gI,o(0)uI,0(0)  =  Cf(0)  +  {x)  uJ'°{x)dx,  I  =  1, . . .  N. 

j= i  Jo 

It  is  worth  pointing  out  that  the  following  result  demonstrates  the  remarkable  difference 
between  nonlocal  quasilinear  hyperbolic  initial-boundary  value  problems  similar  to  (1)  and 
local  quasilinear  hyperbolic  problems  since  it  is  well  known  that  solutions  of  local  problems 
can  attain  discontinuities  in  finite  time. 

Theorem  9  Under  the  additional  assumption  (H9)  the  bounded  variation  weak  solution  of 
problem  (1)  is  continuous. 

Proof.  Let  BT  and  rnr  be  nonnegative  continuously  differentiable  functions.  Furthermore, 
assume  that  g1  is  twice  continuously  differentiable  in  x  and  continuously  differentiable  in 
t.  >  0,  x  G  [0,a’max)  and  gT(xm axp)  =  0,  t  €  [0,  X].  Consider  the  following  initial- 

boundary  value  problem: 

v{  +  {gI{x,t)vI)x  +  mI{xJ.)vI  =  0,  (x,t)  £  (0,xnmx)  x  (0,T), 

^  PXmax  ^ 

<  gf(0,  t)vT(0,  t)  =  CT(t)  +  ^  /  7 1,J/3J(x,t)vJ(x,t)dx,  t  G  (0,T),  (13) 

j=iJo 

,  vT{x,0 )  =  uIfi(x),  x  £  [0, :rmax]. 
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It  follows  from  the  results  in  Section  2  that  problem  (13)  has  a  unique  bounded  variation 
weak  solution.  Furthermore,  using  standard  arguments  (e.g.,  see  [1,  8,  12,  19])  one  can  see 
that  a  characteristic  curve  passing  through  (xj)  is  given  by  (XT(t:  x,t),t),  where  X1  satisfies 

^-x  T{t-xS)  =gI(x :(t;x,t),t) 

and  A"7(t;.x,  t)  =  x.  By  the  assumptions  on  g1  the  function  X1  is  a  strictly  increasing- 
function,  and  therefore  a  unique  inverse  function  tt(x ;  x't)  exists.  Hence  if  we  define  G1  (x)  = 
t1  (x\  0,0),  then  (x,GT(x))  represents  the  characteristic  curve  passing  through  (0,0)  and  this 
curve  divides  the  (x,t)-plane  into  two  parts.  And  the  weak  solution  v1  of  problem  (13)  has 
the  following  implicit  representation 
'  ""’(-Y'fO:  ,'./))• 

exp  {“  Jo  l9x(Xl(s'>  U  *)’  s )  +  m^X^s;  x,  t),  s)]ds}  f  ~  X  ’ 

v^xj)  =  < 

i?7(r7(  0;x,t))- 

exp  {—  f^j(0.x  t)  [gi{XI(s]  x,  t),  s )  +  mI(XI(s]  x,  t),  s)]ds}  t  >  X  ’ 

(14) 


where  RT(t)  =  l/gT(0,t)  C1  (t)  +  Ylj= i  fo™*  7/,J/5J(^,  t)vJ(x,  t)clx  . 

Now,  let  P{t)  =  Ej=1  /0Cmax  w1  (x)  uT(x,t)dx  where  u1  is  the  unique  bounded  variation 
weak  solution  of  problem  (1).  Using  (H2)-(H4)  we  see  that  ^(xp)  =  g:(x,  P(t)),  ^(xp)  = 
PT(x,  P(t ))  and  m1  (x,  t.)  =  m1  (x,  P(t ))  satisfy  the  above  requirements.  Then,  using  (H8),  the 
solution  representation  (14)  and  standard  arguments  as  in  [1]  it  follows  that  v1  is  continuous 
for  this  choice  of  parameters.  Furthermore,  by  uniqueness  of  solutions,  v1  coincides  with 
the  solution  component  u 1  of  the  nonlinear  problem  (1),  /  =  1, . . . ,  N.  This  establishes  the 
result. □ 
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5.  Concluding  Remarks 


In  this  paper  we  presented  a  model  that  describes  the  evolution  of  N  subpopulations 
competing  for  common  resources.  Our  numerical  results  indicate  that  the  parameters  7/,J 
play  a  crucial  role  in  the  dynamics  of  these  subpopulations.  Several  interesting  questions 
about  the  model  (1)  arise  naturally:  What  is  a  good  measure  (in  terms  of  the  rates  g1, 
m1  and  0T)  that  will  lead  to  the  survival  of  the  fittest  in  a  closed  reproduction  case?  In 
an  open  reproduction  case,  which  populations  will  survive  and  which  will  go  to  extinction? 
We  mention  that  for  special  cases  of  structured  age  and  age-size  population  models,  it  was 
proved  in  [17]  that  under  closed  reproduction  a  good  measure  of  species  fitness  is  the  product 
of  the  birth  rate  function  and  the  survivorship  function.  To  our  knowledge,  however,  no 
results  concerning  the  open  reproduction  case  for  structured  populations  are  available.  For  a 
classical  Lotka-Volterra  competition  model  which  is  represented  by  a  system  of  N  differential 
equations,  conditions  on  the  growth  and  mortality  rates  of  each  population  that  will  result 
in  its  survival  or  extinction  have  been  discussed  recently  by  several  researchers  (e.g.,  see 
[5,  6,  7,  21,  22]).  Our  future  efforts  will  focus  on  generalizing  such  results  to  the  distributed 
rate  structured  model  presented  in  (1). 
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